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We present a semi-analytical formulation for calculating the supermodes and corresponding Bloch 
factors of light in hexagonal lattice photonic crystal waveguide arrays. We then use this formulation 
to easily calculate dispersion curves and predict propagation in systems too large to calculate using 
standard numerical methods. 



I. INTRODUCTION 

Coupled optical waveguides are an important class of 
devices for the manipulation and processing of light and 
include, for example, the widely used directional coupler 
[1]. In waveguide arrays, which can be seen as gener- 
alizations of the directional coupler, the propagation of 
light is more complicated. The spreading out of light 
into neighboring waveguides in an array can be seen as 
a discrete version of diffraction. By varying the angle of 
the incoming light, Eisenberg et. al. showed that the 
magnitude and sign of this discrete diffraction could be 
manipulated in ways not possible in regular diffraction 
[2|. The modeling of directional couplers and complex 
waveguide arrays can be achieved with fully numerical 
methods, but the ability to predict unique phenomena 
such as discrete diffraction, is possible only by harness- 
ing the efficiency and physical insight gained by modal 
methods. Although an exact modal method exists, Yeh 
showed that much analytical insight can be achieved in 
the tight binding limit jBj. In this limit the overarching 
coupled waveguide modes or supermodes of the system 
can be represented as linear combinations of the single 
waveguide modes with coupling between adjacent guides 
determined by the overlap between the modal fields and 
the neighboring waveguides. This method shows that for 
uniform waveguides the even superposition is always the 
fundamental coupled waveguide mode. This is consistent 
with the oscillation theorem, according to which, in a ID 
geometry, the mode with the fewest number of nodes has 
the lowest energy [5]. Yeh's method allows us to define 
the symmetry of the individual supermodes and accu- 
rately predicts the evolution of light in the individual 
waveguides within the array. 

The coupling of photonic crystal waveguides (PCWs) 
has been a prominent area of research because of their 
short coupling lengths and unique dispersion properties, 
making them suitable for use in ultra compact devices. 
In square lattice PCW arrays, Locatelli et. al. showed 
that, depending on the number of rows between the wave- 
guides, the coupling coefficient can change sign, inverting 
the sign of the discrete diffraction coefficient throughout 
the array [SI . The ability to invert the sign of the diffrac- 



tion has applications in fields of negative index media 
and imaging. Locatelli's work arises from the work by 
de Sterke et. al. into the coupled modes of square lat- 
tice PCWs [S]. They found that, as the number of rows 
between the waveguides varies, the symmetry of the fun- 
damental mode alternates between even and odd because 
the underlying Bloch mode undergoes a 7r-phase change 
every transverse period. Therefore the order of the modes 
can invert, while still conforming to the oscillation the- 
orem. Using a novel method combining the scattering 
matrix and tight binding methods, Botten et. al. gener- 
alized this formulation to arrays with arbitrary numbers 
of waveguides 0. They found that, as the number of 
rows between the waveguides is increased, the order of the 
supermodes completely inverts, which implies that the 
diffraction coefficient has reversed. In that work a single 
propagating plane wave approximation is made, mapping 
the problem onto that of the homogeneous waveguide ar- 
ray. Through node counting it is shown that the order 
of supermodes in the square lattice PCW array upholds 
the oscillation theorem. 

Despite the significant advantages of hexagonal lattice 
PCW systems, leading to their prominence in experimen- 
tal research, modal methods for hexagonal lattice PCW 
systems are lagging behind their square lattice counter- 
parts. This is because of the analytic complexities which 
arise when modeling the hexagonal lattice symmetry, for 
example, the non orthogonality of the lattice vectors (see 
Fig.Ha)). 

In hexagonal lattice coupled PCWs, the geometry of 
the system can either be in-line or staggered depend- 
ing on the number of rows between the waveguides (see 
Fig. [T]^b)). Ha et al. have shown that this staggered 
geometry is particularly efficient in the coupling of slow 
light between the waveguides [5] . In the staggered geom- 
etry, the coupled PCW modes are degenerate at the Bril- 
louin zone (BZ) edge, where the two approaching modes 
have equal and opposite group velocities. Sukhorukov 
et. al. showed that this is associated with the pres- 
ence of vortex modes, satisfying the Bragg condition f3]. 
In previous work, we show that additional degeneracies 
can occur within the BZ, where the coupled PCW modes 
braid around each other. By analyzing the modes of the 
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cladding, we find a region in the BZ where the two least 
evanescent Bloch modes decay at the same rate, the dif- 
ference in the Bloch mode periodicity leads to the beating 
in field, causing braiding [TU*. 

Generalizing this theory to accommodate an array of 
coupled PCWs is not straightforward. There are a num- 
ber of extra obstacles that arise in modeling this system. 
First, the advantage of using modal method analysis for 
complex problems, such as an array, comes from exploit- 
ing the symmetries of the system, but in the array these 
symmetries are much harder to find. To achieve the re- 
quired symmetry in a system of more than two wave- 
guides, we need to represent the field in each waveguide 
by one quantity which must symmetric about the wave- 
guide center. But in this formulation, we wish to model 
systems where waveguides are created by modifying the 
radius or refractive index of a row of inclusions, rather 
than being completely removed. This means that we can- 
not treat the waveguide as a uniform medium — instead 
we must work completely in Bloch mode bases, in which 
we cannot easily relate the field at the center of the wave- 
guide to its edges. Secondly, when comparing vertically 
aligned waveguide centers in the staggered geometry, we 
find they are at inequivalent points of the unit cell and 
therefore are not related by Bloch's theorem. Thirdly, 
since it has been shown there are two equally important 
Bloch modes in the hexagonal lattice (TU], the treatment 
must be generalized to the vector case, necessitating a 
treatment using matrices rather than scalars. Lastly, the 
basis vectors are not orthogonal which means that the 
Bloch factors in directions parallel and transverse to the 
waveguide do not decouple. 

For this reason we report here a formulation which al- 
lows us to calculate the modes of a arbitrary number of 
coupled PCWs in a hexagonal lattice. We consider the 
tight-binding limit, and, in spite of the complex geome- 
try, find that the resultant matrix has the same symmetry 
as that of Yeh. The key difference is that the coupling 
between adjacent waveguides is expressed in terms of re- 
flection and transmission matrices [7], rather than over- 
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FIG. 1. (a) Unit cell for a hexagonal lattice PC and corre- 
sponding lattice vectors, (b) The two different coupled wave- 
guide geometries in a hexagonal lattice. Staggered (left), sep- 
arated by an even number of defect rows, and in-line (right), 
separated by an odd number of defect rows. 
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FIG. 2. (a) Bloch mode amplitudes defined across an interface 
between PC 1 and PC 2, described in Eq. (b) Bloch mode 
amplitudes defined in a PCW comprised on PC 2 sandwiched 
between layers of PC 1. (c) Shifted Bloch mode amplitudes 
defined in the same PCW as (b). 



lap integrals. This allows for a simple, numerically effi- 
cient, tractable modeling tool which provides the phys- 
ical understanding to exploit the properties of coupled 
hexagonal PCWs. We show, for example, how the braid- 
ing behavior, previously shown for two coupled PCWs, 
generalizes to the PCW array. 

The outline of this paper is as follows: Section [ll] 
presents a detailed derivation of the formulation, showing 
the main result of this paper, before Section |III| shows 
results achieved using the new formulation, comparing 
them to established methods and extending it to con- 
sider problems too complex to be practically calculated 
using standard numerical methods. Section |IV| sums up 
the results of this paper. 



II. MODAL FORMULATION 

In this section we derive a formulation for calculat- 
ing the supermodes of a hexagonal lattice array with an 
arbitrary number of PC waveguides. We begin by deriv- 
ing a resonance condition that may be used to find the 
dispersion relation and modes of a single PCW. We then 
calculate the reflection and transmission properties of the 
barriers between waveguides in the PCW array. Finally 
we use these quantities, together with a tight binding ap- 
proximation, to relate the fields in each waveguide and 
derive a resonance condition that may be used to find the 
dispersion relation and supermodes of the entire PCW 
array. 

We work exclusively in the Bloch mode bases of the 
two PCs that comprise the system of coupled waveguides; 
the Bloch mode basis is the natural basis for a periodic 
medium. At a single frequency, the field in a PC may 
be written as a superposition of propagating and evanes- 
cent Bloch modes [TTl, therefore we represent field by 
vectors of Bloch mode amplitudes. In any given unit 
cell we may write a vector c~ of downward propagat- 
ing/decaying mode amplitudes, and a vector c+ of up- 
ward propagating/decaying modes. 

In each bulk PC, the Bloch amplitudes in different unit 
cells are related by Bloch's theorem. In our notation, the 
first lattice vector ai = {d^ , 0) is parallel to the wave- 
guide. The other lattice vector is a2 = {—dx/2, —dy) (see 
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Fig. [T|a)). Bloch's theorem may be used to define the 
Bloch factor, = e^^''^'^, which is the ratio of the am- 
phtudes of mode s at two points separated by the lattice 
vector afl. Here ks is the Bloch mode's wavevector; in 
our method, its x component is specified for a partic- 
ular frequency, and its y component is found numerically 
using a Rayleigh multipole method 12J. 

Due to an impedance mismatch at the PC's edges, the 
Bloch modes may couple at an interface between PCs 
[13] . Therefore in our calculations we need to consider a 
vector of several Bloch modes in each PC. The number of 
modes required for an accurate calculation varies and is 
related to the number of propagating grating diffraction 
orders; for the examples studied here, two Bloch modes in 
each direction are sufficient, but we include five modes for 
additional accuracy. The unique mode-braiding behav- 
ior of hexagonal lattice PCWs arises from the existence 
of equally dominant evanescent Bloch modes |10J, and 
therefore fundamentally requires multiple Bloch modes 
to be modeled. 

We can consider a PC as a stack of gratings and apply 
the grating diffraction equation. At an interface between 
PCs where the periodicity in x is conserved, coupling 
can only occur between modes with the same This 
means that for an entire system of modes, fiis = 6**^^^=", 
which represents translation by ai, is fixed. On the other 
hand, /i2s is different for each mode. Given a vector 
of downward propagating/decaying mode amplitudes in 
a unit cell, the amplitudes in the unit cell positioned by 
a2 with respect to the first is Ac~, where A = diag(/i2s). 
For a vector c+ of upward propagating/decaying modes, 
the corresponding vector in the unit cell displaced by a2 
is A'c+, where A' — diag(/i2s')- It may be shown using 
reciprocity that for hexagonal lattice PC A'^^ — Ae''^="''^ 
[TT] . When we relate Bloch mode amplitude vectors in 
different unit cells, we generally do so for cells displaced 
by integer multiples of a2. 

Each PC has its own set of Bloch modes. Across an 
interface between PCs (Fig. [2]ja)), the Bloch modes of 
the two PCs couple. Analogously to thin film optics, the 
outgoing Bloch mode amplitudes and are related 
to the incoming Bloch mode amplitudes cj~ and by 
reflection and transmission matrices, which may be cal- 
culated most conveniently from numerically calculated 
PC impedances |13) : 



deriving a result we use in Sec. II C to find PCW array 



C+ ==Ri2Cj- +T21C+ 
C2 = R2lC^ + Ti2Cj^ 



(la) 
(lb) 



supermodes in which the field in the individual PCWs is 
approximately even. The PC waveguide consists of a cen- 
tral waveguiding region (PC 2) sandwiched between two 
PC mirrors (PC 1). The resonance condition is derived 
as a matrix analog to that of an optical slab waveguide. 
We do not assume that the waveguiding region is a uni- 
form dielectric: we allow the possibility that it consists 
of a row of holes — this allows an extra degree of freedom 
when tailoring the dispersion of the PCW. 

We write down the relationship between the Bloch 
mode amplitudes on each side of the waveguide: 



C+ = R21A2C 



R2lA2C^ 



(2a) 



(2b) 



where is the vector of downward mode amplitudes at 
the top of the waveguide, and c"*" is the vector of up- 
ward mode amplitudes at the bottom of the waveguide 
(Fig. ^p). From the up-down symmetry surrounding the 



interface, we know = R2i- 



The asymmetric factor e*''^''^ arises due to the non- 
orthogonality of the two lattice vectors. Its presence be- 
comes particularly inconvenient in Sec. |II C[ which re- 
quires a high degree of symmetry. We remove this factor 
by defining 



(3) 



where xrow = {vr/ dy){dx/2). The term xrow corre- 
sponds to the shift in x associated with a translation by 
{—VrI dy)ai2, where ur is the y coordinate of the unit cell 
for which is defined. Note that the c+ and c~ in 
Eqs. ([2]) are defined on opposite sides of the waveguide 
region (Fig. [2jb)), and so their xrow differ by d^/^: this 
is what removes the symmetry-breaking quantity e^^^'^^ , 
The original quantities are defined in unit cells sepa- 
rated by multiples of the diagonal lattice vector a2. For 
every second row of unit cells, multiplication by e~*'^^^ 
corresponds to a translation by an integer multiple of ai , 
and so has the straightforward interpretation of be- 
ing the amplitude vector of the Bloch modes in the unit 
cell with a; = in that row. For the other rows, has 
no direct interpretation, but c^Q^k:,d^/2 -j-j^g vector of 
Bloch mode amplitudes at x = dx/2. This new basis for 
Bloch mode amplitudes necessitates the definition 



Ae' 



(4) 



A. The single waveguide problem 

First, we derive the resonance condition for a single 
PC waveguide, which may be used to find the PCW's 
dispersion relation and modes. This provides an intro- 
duction to our nomenclature and methods, as well as 



where A replaces A in relating the downwards amplitudes 
of Bloch modes in different rows of a PC. The reflection 
and transmission matrices are unchanged. 

With our new quantities, we can rewrite Eqs. ^ as 



c+ = R21A2C , 

= R21A2C+, 



(5a) 
(5b) 



relating shifted Bloch mode amplitudes (Fig. ^c)) By 
substituting Eq. ( 5b ) into Eq. ( 5a I , we find the resonance 



condition for a single PCW mode, 

= R2lA2R2lA2C^, (6) 

0- (I±R2iA2)c-. (7) 

So, the condition for a single waveguide mode is 

det(I±R2iA2) = 0. (8) 

The eigenvector c~ associated with the zero eigenvalue 
gives the waveguide mode in terms of the Bloch modes of 
PC 2, the in-band PC. Using Eqs. ([5]) and ([T]), we derive 



c"*" = +c , or 
c+ = c , 



(9a) 
(9b) 



where Eq. ( |9a| ) holds for even PCW modes and Eq. ( |9b[ ) 
holds for odd PCW modes. In an experiment it is easier 
to couple light into symmetric modes, so we will concen- 
trate on the even PCW modes. 

The contribution of even modes to the field in a wave- 
guide may be summarized as 



1 



(10) 



which may be shown to represent the amplitudes of the 
even superpositions of Bloch modes in a single waveguide 
(Eq. ([12|). 

To find the waveguide mode's field H(x, y) inside a 
waveguide unit cell, we convert our calculated result, c+ 
and c~ , back into the basis of c+ and . If we choose c~ 
such that it is defined with respect to the origin, then by 
Eq. ^ there is no difference between the two bases, i.e., 
= c~ and the amplitudes of the upward Bloch modes 
in this cell are A2C+. Therefore the magnetic field Hs^g 
at position r in the unit cell is 



i7swg(r) = (A2C+)^*+(r) +c-^*-(r), 
= e^(A2*+(r) + *-(r)) 



(11) 
(12) 



where ^^^(r) and ^E'~(r) are respectively vectors of the 
fields of the forward and backward Bloch modes. A sim- 
ilar expression may be written for the electric field. It 
is straightforward to use the transmission matrix T21 to 
calculate the Bloch mode amplitudes of PC 1 from c+, 
and a similar procedure may be followed to calculate the 
field of the waveguide mode throughout PC 1. 



B. The transmission and reflection matrices across 
the waveguide barrier 

In this paper we consider PCW arrays consisting of 
layers of two different photonic crystals, PC 1 and PC 2, 
where PC 1 is the bulk PC, which is in bandgap, and 
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FIG. 3. (a) Bloch mode amplitudes defined across a barrier 
consisting of £ rows of PC 1 between two waveguides. The 
broken lines indicate a perfectly matched layer, rather than 
an interface. (b) Bloch mode amplitudes defined in a section 
of the PCW array. 



PC 2 is the waveguide . The barrier consists of £ rows 
of PC 1, and the waveguide is one row thick. PC 1 and 
PC 2 may be dielectrics or PCs, and need not have the 
same background refractive index. 

This section focuses on a part of the array consist- 
ing of just two neighboring waveguides and the barrier 
separating them (see Fig. |3^), to find the reflection and 
transmission matrices across the barrier. We neglect re- 
flections off the outermost edge of each waveguide, and 
aim to relate the vectors cj, c^, cj^-^ and c'j'^^ to find 
the reflection and transmission matrices in each direction 
across the barrier, 

c+ = RA2CT + T'A2C+_i, (13a) 

cT+i - TA2C7 + R'A2C++i, (13b) 

where R, T and R',T' are reflection and transmission 
matrices from above and below the barrier respectively. 

If we set c^_|_j = (see Fig. [3^), then the only Bloch 
mode incident on the barrier is from above, i.e. cJ, and 
our relations in Eqs. (13a I and (13b) become = -^^^7 



and Cj^^ = Tc^ . The relations between and c^^^ 
and the Bloch modes inside the barrier, are given by 



c+ =R2iA2C- +T'i2A^iC+, 



c+ 



Ti2A^iC,-, 

Rl2A^lC^, 

■R' \^r+ - 



T21A2C+ 1. 



(14a) 
(14b) 
(14c) 
(14d) 



From the up-down symmetry, we know T'j^j — T12 and 



R'l2 



R12. We can rearrange these relations to find 



Cj_^_i in terms of A2C^- : 

=Ti2Al(I- Ri2A^iRi2A0-'T2iA2C- = TA2C-, 

(15) 

^ T =Ti2A((I - Ri2A^iRi2Af )-iT2i. (16) 
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Similarly, by rearranging the relations to find in terms 
of A2C~, we derive 

R = R21 + T12 A^iRi2 A^i (I- R12 R12 A^i )-iT2i . (17) 



If we instead set cj = 0, then Eqs. (13a) and (13b) 
become = T'c^_^j^ and cj_^_-^ = R'c^_)_j^. By following 

the same method as above, we derive T" = T and R' — 
R. 



C. PCW array dispersion relation and supermodes 

In a waveguide array, interactions occur between ev- 
ery pair of waveguides. For large arrays of many wave- 
guides, it is unfeasible to model all these interactions. 
Following Yeh we ignore the interactions between 
non-adjacent waveguides, making a nearest neighbor ap- 
proximation, which is valid in the tight binding regime. 
In doing so, we assume that the field decays substantially 
across the barriers between waveguides, i.e., we assume 
that l/i^l <C 1 for every element fi of Ai. Formally, we 
neglect terms of order 0(A^^). Doing so excludes the 
interactions between non-adjacent waveguides, as well as 
the Fabry-Perot-like terms in the interactions between 
adjacent waveguides. 

Under this approximation, the transmission and reflec- 
tion matrices in Eqs. (16) and (17) simplify to 



T = Ti2AlT2i, 



R — R21- 



The reflection matrix for the waveguide barrier is ap- 
proximated by the reflection off an infinitely thick barrier 
since any effect from the neighboring waveguide consists 
of terms of order O(Af^). 

We can now relate the Bloch mode amplitudes of each 
of the waveguides in an Af-waveguide PCW array. The 
vectors c^'" and cJ that represent upward and downward 
Bloch mode amplitudes in waveguide j are as defined 
in the previous section. The topmost and bottommost 
waveguides, 1 and M, interact with only one neighbor 
each. The relationship between and c'^ is similar to 
that of a standalone PCW (Eq. ([5])); the difference is that 
Cj^ includes a contribution TA2CJ from waveguide 2: 



= R2lA2C]^, 

c+ = R21A2CJ; + TA2C, 



(19a) 
(19b) 



Similarly, in waveguide Af , c^j 
from waveguide M — 1: 

Cm = R21A2C+ + 
Cm = R2iA2Cjy^. 



includes a contribution 



(20a) 
(20b) 



The other waveguides 1 < j < Af, interact with both 
their neighbors (see Fig. ^ 



(21a) 
(21b) 



c^- = R21A2C+ + TA2C-_i, 

C+=R2lA2C-+fA2C++l. 



In these equations, we represent the field in each wave- 
guide by two different quantities: on the top edge as the 
amplitude vector of the forward modes c~ , then on the 
bottom edge as the backward modes c+. In order to 
apply methods analogous to those used with coupled di- 
electric waveguides [3], we need to represent the field in 
each waveguide by one quantity. We define 



(22) 



Gj only represents the contributions of even superpo- 
sitions of Bloch modes in the waveguide, as noted in 
Sec. |II A[ From here on we represent the amplitudes in 
waveguide j by . We start by using Eq. ( 22 ) to write 



the 2M equations (19)-(21) as A/ equations relating the 



ei = RA2ei + ^TA2C+, 
e, = RAse, + ^tA2C++i + ^TA2C-_i, 
em = RA2eA/ + ^'TA2Cm_-^. 



(23a) 
(23b) 
(23c) 



(18a) order to replace the residual in the last terms of 

these equations, we recall our approximation discard- 
ing terms of order 0(A^^). The field in each wave- 

(18b) 

guide is an even single waveguide mode, plus interac- 
tions of order O(A^) from adjacent waveguides: thus 

we write 



for each waveguide in the array, using Eq. ( 9a 



c+ = Cj +0{A.{). Note that the cf terms in Eq. ^ are 
all multiplied by T, which is of order 0{A{). Therefore 
for these terms we can ignore the 0{A.{) contributions 



from other waveguides and write c 

Applying this substitution to Eq. (23), we write the 
equations solely in terms of e-,-, in the form of a tri- 
diagonal block matrix equation. 



/ABO 
BAB 
B A 



\0 



■ ° \ 









62 


62 







= 63 


• A / 


VeAf/ 


VeA// 



(24) 



where A = R21A2 and B = ^TA2. 

This is the resonance condition for a coupled waveguide 
array, in the same way that Eq. ^ is the resonance con- 
dition for a single waveguide. For each frequency, a range 
of kx can be scanned for the existence of supermodes by 
testing whether the condition may be satisfied by each 
individual k^.. 
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The difference between Eq. (24 1, and the expression 
for conventional waveguide arrays [3,, is that each ele- 
ment in Eq. (24) is a n x n matrix rather than a scalar. 



If n Bloch modes are considered in each waveguide, then 
the matrix in Eq. (24) is of dimension nM. For a large 



number of waveguides M, it can be computationally ex- 
pensive to solve the eigensystem of this matrix. There- 
fore, in a block-matrix analog to the conventional treat- 
ment [3], we exploit the tridiagonal block symmetry of 
the matrix and block-diagonalize it analytically using a 
similarity transformation. The problem then reduces to 
solving the eigensystems of M matrices, each with dimen- 
sion n. Analogous to Yeh's treatment [3], the elements of 
the symmetric and unitary similarity transform $ that 
block-diagonalizes the tridiagonal matrix are given by 



where I has dimension n and 



Af + 1 



M +1 



(25) 



(26) 



To demonstrate how to use Eq. ( 24 ) to calculate a dis- 



persion relation, we give the example of the two wave- 
guide case, for which 



A B 
B A 



(27) 



This can be block diagonalized using a similarity trans- 
formation to give 



with all other Xj = 0, where as = 2cos(s7r/(M -|- 1)). 
[3] . The fields in each waveguide for supermode s can be 
found by 



$ t;X t; 



(31) 



where denotes block column s of $ 

The procedure to find the dispersion relation of a cou- 
pled PCW array therefore involves calculating A and 
B for frequency-fca, pair and checking whether det(I — 
(A — (TsB)) = for any ag. At a frequency and kj. for 
which this determinant is zero, then the null vector of 
I — (A — (JsB) is the vector x^. 

The resonance condition Eq. ( 30 ) can be written di- 



rectly in terms of the original reflection, transmission and 
propagation matrices, using Eq. (|4]): 

(I-e*'==='^/2j^2iA2/+ye^'=^'^(^+i)/2T^2A^i/T2iA2/)x, - 0. 

^ . (^^^ 
The first two terms in Eq. ( 32 ) give the resonance con- 



dition of the single waveguide in a hexagonal lattice. 
Therefore the third term, of order A{, is the effect of 
the neighboring waveguides in the tight-binding approxi- 
mation. In the photonic band gap, all Bloch modes in the 
cladding are evanescent and have Bloch factors < 1. 
As we increase the distance between waveguides, i in- 



creases, and A^ 



0. Therefore the resonance condition 



of the array approaches that of the single waveguide as 
we increase the distance between the waveguides, as ex- 
pected. 



A B 
B A 



(28) 



D. PCW array modes and propagation 



where, from Eq. (26) ^' 



1 

V2 



J ii ) , giving 



A + B \ pi +e2\ /ei +e2\ /ki 
A - By l^ei - 62/ \e1-e2J ~ 1^x2^ 

(29) 

Equation ( |29[ ) decouples into two nontrivial solutions: 
(I — (A -f- Bjjxi = with X2 = (even supermode), 
and (I — (A — B))x2 = with xi = (odd supermode). 
The coupled waveguide array's dispersion relation may 
be constructed from the frequency-fca, pairs at which at 
least one of these conditions is satisfied. So the condi- 
tion for an even supermode is det(I — (A + B)) = 0, and 
for this mode ei = §2 = 5X1, where xi is an eigenvec- 
tor of I — (A + B) with an eigenvalue of zero. This is 
the Bloch mode analogue of the result achieved in our 
previous work 10' . 

Due to the block tridiagonal symmetry of the matrix 
in Eq. ( 24 ) , similar conditions may be obtained even for 
very large systems of coupled waveguides. Generally, the 
condition for the sth PCW array supermode is 



At the kx associated with supermode s, Eq. (30) may 
be used to calculate Xg. The other x^ = 0, and these 
values are used implicitly in Eq. (31) to find the super- 
mode's field amplitudes in each waveguide j. 

In choosing to represent the field in each waveguide by 
the vector , we assumed that the field is an even super- 
position of upward and downward Bloch modes in each 
waveguide. In Yeh's treatment of coupled waveguides [3], 
the field in the waveguide array is written as a superpo- 
sition of single waveguide modes, with one amplitude per 
waveguide. We do the same by writing Xs ~ e, where e 
is the vector of Bloch mode amplitudes calculated for a 
single waveguide. This approximation is consistent with 
an eigenvalue perturbation theory approach: k^ is found 
including terms of order O(A^), but the modal field is 
found ignoring these terms. 

From Eq. (31 ), it follows that the field of supermode s 



in waveguide j is 



(33) 



(A - cTcBlx. 



(30) 



We now generalize to consider a superposition of super- 
modes with amplitudes Ws, and calculate the resulting 



7 



field in each waveguide, ej = a^e, which we simply rep- 
resent by the scalar amplitude dj. Then 



dj ^y^^(f)jsWs, 

s 

i.e. a =(/)w, 



(34a) 
(34b) 



where a and w are vectors of dj and respectively. 

Propagation in the PCW array is most simply rep- 
resented in the basis of its supermodes, each of which 
simply acquires phase as it propagates. We use the 
propagation constant of each supermode s to write 
its Bloch factor /ii^ — exp(ik,j;dx) Given an initial field 
w(0) — (f>~^a{0) at the start of the waveguide array, the 
amplitude in waveguide j after propagating p periods is 



s 

i.e. a{p) =0LPw(O), 



(35a) 
(35b) 



where L = diag(/iis). Writing the initial condition w(0) 
in terms of waveguide rather than supermode amplitudes, 
we see that the field at any integer number of periods p 
along the PCW array is 



sl{p) = (j)L (/}-'- sl{0) 



(36) 



where 4> = 4>~^. 

Therefore, knowing only the dispersion relation, we can 
analytically calculate the amplitudes a(p) of each wave- 
guide mode at an arbitrary number of periods p into the 
structure. The actual field in each waveguide may be 
calculated from a(p), and the field iJswg(r) of the single 
waveguide mode over the domain < x < d. Recall that 
m Eq. (§, we transformed the Bloch mode amplitudes c 
into a mathematically convenient basis and wrote them 
as c. Unlike in Sec. |II A[ we cannot set the origin of the 
coordinate system such that aj = dj for all waveguides 
j, therefore we must explicitly perform a transformation 
into a more natural basis. We map each dj to the ampli- 
tude Uj of the single mode waveguide in the unit cell p 
periods into the structure, by writing 



(37) 



where Kj — exp{ikx(xR — pd)), with xn the coordinate of 
the unit cell in waveguide j (more specifically, the point 
in that unit cell marked with a cross in Fig. [l|^a) . p pe- 
riods into the structure. The origin of coordinates may 
be chosen such that xr = pd for every unit cell p along 
waveguide j — 1: then Kj = 1 for all waveguides that 
are in line with waveguide 1, and Kj — exp{ikxd/2) for 
waveguides staggered by Ax = d/2 with respect to wave- 
guide 1. Here, kx is approximated to that for the single 
waveguide mode. 

Finally, we write the field in the waveguide array as 
a superposition of translated single waveguide modes 
-ffswg(r ~ R-j) and their amplitudes a(p): 



H{v) 



j 



KjdjH^^g(^v 



(38) 



where Rj = (xR^yn) is the coordinates of the nearest 
unit cell to r inside waveguide j. 

If we neglect the evanescent tails of single waveguide 
modes, then we can write the field in or near waveguide 
j of the array in terms of one waveguide mode only: 



Hj{r) = KjdjHs„g{r -Rj). 



(39) 



If we confine our interest to the field inside the wave- 
guide at a single point per unit cell, then r — Rj is con- 
stant and only one value i?swg(r ~ Rj) is used; the mode 
may be renormalized such that iJswg(r — Rj) — 1- Then 
the field throughout the waveguide array is given by the 
elements aj of a(p) , which can be calculated analytically. 



III. RESULTS 

We now apply the methods we have developed to calcu- 
late dispersion relations for arrays of coupled waveguides, 
as well as the propagation of light through them. Where 
the use of conventional numerical tools is practical, we 
compare the results of our method to these. For simplic- 
ity, in all our examples we consider arrays constructed 
from the same two PCs. Both PCs are regular hexagonal 
lattices of cylindrical inclusions with radius r = 0.3 d, 
in a dielectric background with n = 3. The difference 
between the PCs is that the barrier, PC 1, has air holes 
with n — 1, whereas the waveguide material, PC 2, has 
inclusions with n = 1.5. Our method also applies to 
more conventional Wl PC waveguides, and can easily be 
extended to dispersion engineered waveguides: we con- 
sider an example with infiltrated holes in PC 2 in order 
to demonstrate that our method works even when the 
waveguiding region is not a uniform dielectric. 

We consider two arrays based on this waveguide. The 
first is a system of Af = 4 waveguides, each separated by 
a barrier oi £ = 4 rows of PC 1. The second array is a 
system of M — 31 waveguides, also separated by £ = A 
rows of PC 1. The second array involves a structure with 
a width of over 150 unit cells, which is not feasible to sim- 
ulate using conventional numerical tools. The separation 
oi £ — 4 means that the waveguides are staggered. 

We consider 5 Bloch modes in each PC, which means 
that for each frequency-Zca, pair, testing for the existence 
of a supermode only involves calculating the determi- 
nants of M = 4 different 5x5 matrices (Eq. ([30])). At 
the frequencies we consider, we could get results of ac- 
ceptable accuracy using only 2 Bloch modes. 



In Sec. Ill A we use Eq. (30 1 to calculate the disper- 
sion relations of these structures. Then, in Sec. |IIIB] we 
calculate the discrete diffraction pattern of the field as it 
propagates through the waveguide array. 



A. PCW Array Dispersion relations 

First, we calculate the dispersion curves of the super- 
modes of the 4 waveguide system described above, for 
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light polarized with the H field out of the plane, at 
normalized frequencies d/X £ (0.29,0.30). The disper- 
sion curves are shown in Fig. |4] At most frequencies 
there are four supermodes based on the even single wave- 
guide mode, and their dispersion curves agree well with 
those calculated using the Fictitious Source Superposi- 
tion (FSS) method, which is considered to be highly ac- 
curate for this kind of problem [T?] . 




2.4 2.6 2.8 3.0 

kxd 

FIG. 4. Comparison of the dispersion relation calculated us- 
ing the FSS (green dots) , and our semi-analytic modal method 
(red lines) for a four waveguide array with waveguides sepa- 
rated by four rows of inclusions (shown in inset). Background 
refractive index, n = 3, cylindrical inclusions of index n = 1 
with the waveguide created by changing the refractive index 
of the waveguide inclusions to n = 1.5. 

Having established the accuracy of our method with 
a relatively small array of 4 waveguides, we now show 
its power by finding the dispersion curves of the super- 
modes of a 31-waveguide array (Fig. [5]). We see that 
for this increased number of waveguides and supermodes, 
the dispersion curves of the supermodes start to form a 
band, but the dispersion curve still exhibits the braiding 
behavior shown for the coupled PCW system pD]. 




0.290 



2.0 2.2 2.4 2.6 2.8 3.0 

kxd 

FIG. 5. The dispersion curve for a M = 31 PCW array with 
four rows of inclusions between each waveguide. The back- 
ground index is n = 3, the bulk photonic crystal inclusions 
have index n = 1, and the waveguide is formed with inclusions 
of index n = 1.5. 



B. Driven PCW Arrays 

We now turn our attention to field propagation 
through the two PCW arrays. We seek to simulate a 
waveguide array where the field is incident from a sin- 
gle waveguide. In the M = 4 waveguide array, we set 
a(0) = (0,0,1,0)"^, corresponding to an initial condition 
where the entire field is in the 3rd waveguide. At a nor- 
malized frequency of d/A = 0.293, there are four prop- 
agating supermodes, at k^d = —2.68, —2.70, —2.74, and 
—2.80. A single waveguide with these parameters has 
kxd = —2.72 We choose the negative values of kx as they 
are associated with positive group velocity (Fig. |4]). We 
calculate their field at one point per unit cell over the 
first 100 periods of the waveguide array, i.e., we calculate 
a(p) from p = to p — 100 (Fig. [6]). The approximation 
made when writing ~ e has normalized residual errors 
of (2.2,0.3,0.3,2.1) X lO'^ for the four supermodes. 

These results are indistinguishable from the results of 
an FSS simulation, although neither method simulates 
the coupling at the interface between the input waveguide 
and the waveguide array. In order to include coupling 
effects, we use a well established Bloch modal supercell 
scattering method [IV to model propagation through a 
single waveguide into the waveguide array. We sample 
the field at the same points in each unit cell for each 
waveguide in the array. There is good agreement between 
the supercell method and our results (Fig. |6]). 




20 40 . 60 80 100 

Period 



FIG. 6. A comparison of methods for calculating the field 
throughout the waveguide array at frequency d/X = 0.293. 
One point is taken per unit cell of the array shown in Fig. [4] 
The colors represent the different waveguides; the open mark- 
ers are from Bloch mode supercell calculations (SCM) and the 
solid markers are results from our modal method (MM). 

Now we consider the system of M — 31 coupled wave- 
guides. Again, we start with all the field in the cen- 
tral waveguide, so all elements dj of a(0) are 0, ex- 
cept flig — 1. We calculate the propagation of the field 
through 1100 periods of the waveguide array (Fig. [t]). 
For the first 450 periods, a typical discrete diffraction pat- 
tern is observed. After around 500 periods, the diffracted 
light reaches the outer waveguides, and edge effects start 
to play a significant role in the field profile. Once the 
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kx values of the supermodes are calculated, generation of 
the data for Fig. [7] takes under 10 seconds on a desktop 
computer. 



IV. DISCUSSION AND CONCLUSION 

The semi-analytical formulation developed here pro- 
vides an accurate modal method to solve for the modes of 
coupled PCW waveguides. An advantage of this method 
is that the computational time increases linearly as the 
number of waveguides increases, in the sense that for M 
waveguides, M equations of the form of (32 1 need to be 



solved to obtain the PCW modes at each frequency. The 
advantage here is that each mode is identified by its sym- 
metry and solved for independently. Therefore, the mem- 
ory requirements remain constant as M increases. This 
contrasts with supercell methods, where as the number 
of waveguides increases, the computation domain grows, 
increasing both the required memory and computation 
time. Furthermore, in supercell methods the modes of 
different symmetry are not readily identifiable, making 
crossings in the dispersion curves, such as those shown 
in Figures |4] and [5] difficult to recognize. 

Our method is, in essence, a tight-binding formulation 
and therefore it warrants comparison to a conventional 
tight-binding formulation based on overlap integrals [3]. 
In such conventional formalisms one solves for the iso- 
lated single PCW mode and assumes that the modes of 
the coupled PCWs can be written in terms of a super- 
position of isolated single PCW modes. The interaction 
between adjacent PCWs is found by computing overlap 
integrals. The formulation presented here is identical in 
terms of the assumptions it is based on, i.e. adjacent 
waveguides couple weakly, however in place of overlap 
integrals interaction between waveguides is modeled by a 
transmission matrix. There is one key advantage in using 
our modal method based calculation: we find the PCW 
modes by finding the propagation constants kx satisfying 
(32) at a given frequency. In contrast, an overlap inte- 



gral based formulation for PCWs gives the shift in fre- 
quency of supermodes from the single PCW mode. The 
former is much more desirable for modeling propagation 



problems, as once an operation frequency is chosen the 
relevant modes can be obtained in a single set of cal- 
culations. At the level of Bloch modes and supermodes, 
the similarities between the acoustic and electromagnetic 
treatments mean that this modal method may also be ap- 
plied to arrays of elastic waveguides [121 HB] . 

In the work of de Sterke [S], it is shown that al- 
though the fundamental mode can be either an even or 
odd superposition of PCW modes, when taking a cross 
section along the PCWs, the fundamental mode always 
possesses one less node than the second mode. Here we 
find that, when taking a cross section along the coupled 
PCWs, node counting cannot be used to predict the or- 
der of modes in the hexagonal lattice. This is because 
the oscillation theorem only holds for one dimensional 
Sturm-Liouville problems. In the case of square lattice 
PCWs, the modes have a frequency below the first Wood 
anomaly 6 . This allows for the use of a scalar approach 
where coupled PCWs are mapped on to coupled uniform 
waveguides. The rich physics found in hexagonal lattice 
PCW arrays, is due to the fact that there are two domi- 
nant Bloch modes and a vectorial approach is needed. We 
have observed that near the center of the Brillouin zone, 
where one operates below the first Wood anomaly, the 
oscillation theorem holds in the hexagonal lattice. How- 
ever, as the Bloch wavevector approaches the Brillouin 
zone edge, crossing the first Wood anomaly and espe- 
cially upon entering the braided region, the modes be- 
come inherently vectorial in nature, ic. their Bloch mode 
expansion requires at least two Bloch modes, and thus 
cannot be mapped on to a scalar amplitude like in the 
square lattice. Therefore there is no general node count- 
ing algorithm governing the order of modes in hexagonal 
PCW arrays. 
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